clear
version 14.0
set more off
set scheme s2mono
capture log close
log using "..\logs\lasso - select soils.txt", text replace

* Specify the year of cash rent data to use in the regression
local reg_year 2013
* Specify 1 if irrigated or 0 if nonirrigated
local practice 0

use "..\dataAnalysis\rent_climate_data", clear

if `practice'==0{
gen ln_cash_rent=ln(cash_rent_non`reg_year')
}
if `practice'==1{
gen ln_cash_rent=ln(cash_rent_irr`reg_year')
}

*------------------------------------------------------------------------
* Potential soils
*------------------------------------------------------------------------
local soils  slope om sar gypsum cec7 ec  ksat clay_perc silt_perc ///
			bulkDensity soc0_150 rootznemc ph_less6 ph_greater7dot5 hydgrp* text_*

local soils_nonlin  slope_* om_* sar_* gypsum_* cec7_* ec_*  ksat_* clay_perc_* silt_perc_* ///
			bulkDensity_* soc0_150_* rootznemc_* ph_less6_* ph_greater7dot5_* hydgrp* text_*

*------------------------------------------------------------------------
* Controls to include in the regression
*------------------------------------------------------------------------
local deficit water_deficit
local surplus water_surplus
local edd dday34C

local knots_d 2
local knots_s 2
local knots_g 3
local knots_e 3

if `knots_d'==2 {
capture drop `deficit'_spline*
gen `deficit'_spline1 = `deficit'
matrix `deficit'_knots=[0]
}

if `knots_d'>=3 {
capture drop `deficit'_spline*
mkspline `deficit'_spline = `deficit', cubic nknots(`knots_d') 
matrix `deficit'_knots=r(knots)
}

if `knots_s'==2 {
capture drop water_surplus_spline*
gen `surplus'_spline1 = `surplus'
matrix `surplus'_knots=[0]
}

if `knots_s'>=3 {
capture drop water_surplus_spline*
mkspline `surplus'_spline = `surplus', cubic nknots(`knots_s') 
matrix `surplus'_knots=r(knots)
}

if `knots_g'==2 {
capture drop gdd_spline*
gen gdd_spline1 = gdd
matrix gdd_knots=[0]
}

if `knots_g'>=3 {
capture drop gdd_spline*
mkspline gdd_spline = gdd, cubic nknots(`knots_g') 
matrix gdd_knots=r(knots)
}

if `knots_e'==2 {
capture drop `edd'_spline*
gen `edd'_spline1 = `edd'
matrix `edd'_knots=[0]	
}

if `knots_e'>=3 {
capture drop `edd'_spline*
mkspline `edd'_spline = `edd', cubic nknots(`knots_e') 
matrix `edd'_knots=r(knots)		
}	

*******************
* Step 1
*******************

local knots_`deficit'm1=`knots_d'-1
local knots_`surplus'm1=`knots_s'-1
local knots_gddm1=`knots_g' -1
local knots_`edd'm1=`knots_e'-1

tab stfips, gen(state_fe)
* Drop one state to omit as baseline
drop state_fe22

forvalues j=1/`knots_`deficit'm1' {
	lassoShooting `deficit'_spline`j' `soils_nonlin', controls(state_fe*) verbose(0) fdisplay(0)
	local `deficit'`j'Sel `r(selected)' 
	di "``deficit'`j'Sel'"
}

forvalues j=1/`knots_`surplus'm1' {
	lassoShooting `surplus'_spline`j' `soils_nonlin', controls(state_fe*) verbose(0) fdisplay(0)
	local `surplus'`j'Sel `r(selected)' 
	di "``surplus'`j'Sel'"
}

forvalues j=1/`knots_gddm1' {
	lassoShooting gdd_spline`j' `soils_nonlin', controls(state_fe*) verbose(0) fdisplay(0)
	local gdd`j'Sel `r(selected)' 
	di "`gdd`j'Sel'"
}

forvalues j=1/`knots_`edd'm1' {
	lassoShooting `edd'_spline`j' `soils_nonlin', controls(state_fe*) verbose(0) fdisplay(0)
	local `edd'`j'Sel `r(selected)' 
	di "``edd'`j'Sel'"
}

*******************
* Step 2
*******************
lassoShooting ln_cash_rent `soils_nonlin', controls(state_fe*) verbose(0) fdisplay(0) 
local ySel `r(selected)'
display "`ySel'"

*******************
* Step 3
*******************
* Get union of selected instruments
local SelSoils `ySel' 
forvalues j=1/`knots_`deficit'm1' {
	local SelSoils : list SelSoils | `deficit'`j'Sel
}
forvalues j=1/`knots_`surplus'm1' {
	local SelSoils : list SelSoils | `surplus'`j'Sel
}
forvalues j=1/`knots_gddm1' {
	local SelSoils : list SelSoils | gdd`j'Sel
}
forvalues j=1/`knots_`edd'm1' {
	local SelSoils : list SelSoils | `edd'`j'Sel
}
di "`SelSoils'"
global SelSoils `SelSoils'

reg ln_cash_rent `deficit'_spline* `surplus'_spline* gdd_spline* `edd'_spline* $SelSoils state_fe* ///
		[aw=croplandNon_acres], vce(cluster stfips_clust)

* This code saves $SelSoils which can then be used to specify soil controls in my other code
* to estimate damages

log close
